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Abstract 



We describe a Monte Carlo scheme for simulating polydisperse fluids 
within the grand canonical ensemble. Given some polydisperse attribute a, 
the state of the system is described by a density distribution p{a) whose form 
is controlled by the imposed chemical potential distribution /x(cj). We detail 
how histogram extrapolation techniques can be employed to tune /^(ct) such as 
to traverse some particular desired path in the space of p{(7). The method is 
applied in simulations of size-disperse hard spheres with densities distributed 
according to Schulz and log-normal forms. In each case, the equation of state 
is obtained along the dilution line, i.e. the path along which the scale of p{a) 
changes but not its shape. The results are compared with the moment-based 
expressions of Boublik et al (J. Chem. Phys. 54, 1523 (1971)) and Salacuse 
and Stell (J. Chem. Phys. 77, 3714 (1982)). It is found that for high de- 
grees of polydispersity, both expressions fail to give a quantitatively accurate 
description of the equation of state when the overall volume fraction is large. 
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I. INTRODUCTION 



Statistical mechanics was originally formulated to describe the properties of systems of 
identical particles such as atoms or small molecules. However, many materials of industrial 
and commercial importance do not fit neatly into this framework. For example, the particles 
in a colloidal suspension are never precisely identical to one another, but have a range of 
radii (and possibly surface charges, shapes etc). This dependence of the particle properties 
on one or more continuous parameters is known as polydispersity. 

To process a polydisperse colloidal material, one needs to know its phase behavior, i.e. the 
conditions of temperature and pressure under which a given structure is thermodynamically 
stable. The main obstacles to gaining this information arise from the effectively infinite 
number of particle species present in a polydisperse system. Labeling these by the continuous 
polydispersity attribute a, the state of the system must be described by a density distribution 
p((t), rather than a finite number of density variables. The phase diagram is therefore infinite 
dimensional, a feature that poses serious problems to experiment and theory alike. 

The chief difficulty faced in experimental studies of polydisperse systems is that the 
infinite dimensionality of the phase diagram precludes a complete mapping of the phase 
behavior. Instead one is forced to focus attention on particular low dimensional manifolds 
(slices) of the full diagram. Typically this involves determining the system properties along 
some desired trajectory through the space of p(o"). Such a strategy is often pursued in 
experiments on colloidal suspensions where the phase behavior is studied along a so- 
called dilution line. The experimental procedure for tracking this line involves adding a 
prescribed quantity of colloid of some known degree of polydispersity to a vessel of fixed 
volume V, the remaining volume being occupied by a solvent. The number distribution 
of colloidal particles N{a) determines the density distribution of the suspension, p(cr) = 
N{a)/V. Since in a given substance, the relative proportions of the number of particles of 
each a are fixed, changing the amount of colloid added simply alters the scale of p(cr), not 
its shape. Thus, by varying N{a) at fixed V (or vice versa) one traces out a locus in the 
phase diagram in which only the overall scale of p{a) changes. 

As regards theoretical studies of phase behavior, these typically endeavor to calculate the 
system free energy as a function of a set of density variables. The difficulty in achieving 
this for a polydisperse system is that the free energy /[p(cr)] is a functional of p(cr), and 
therefore itself occupies an infinite dimensional space. This renders intractable the task of 
identifying phase boundaries. Recently, however, progress has stemmed from the observation 
that it is possible to approximate the full free energy by a so-called "moment free energy" 
containing the full ideal gas contribution plus an excess part that depends only on certain 
principal moments of the full excess free energy ||^. Doing so reduces the problem to a finite 
number of density variables and allows calculation of phase coexistence properties within 
a systematically refinable approximation scheme. Additionally the theory delivers (for the 
given free energy) exact results for the location of spinodals, critical points and the cloud and 
shadow curves. Use of this approach promises to enhance significantly our understanding of 
phase behavior in polydisperse systems. 

In view of the approximations inherent in theoretical approaches, it is natural to consider 
deploying computer simulation to study the phase behavior of polydisperse colloidal fiuids. 
Although simulations (like experiment) are restricted to studying limited regions of the phase 
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diagram such as the dilution hue, they have the advantage that can be used to investigate 
the same model systems as studied theoretically. Furthermore, they deliver (modulo finite- 
size effects) essentially exact results, providing invaluable benchmarks against which to test 
theoretical predictions. Sometimes too, the physical insight gleaned from simulations serves 
as the impetus for fresh theoretical advances. 

One simulation approach for obtaining the thermodynamic properties of a polydisperse 
system is to simply mimic the experimental procedure. This can be achieved by employing 
a canonical ensemble (CE) formalism wherein a simulation box of fixed volume V is pop- 
ulated by a prescribed number of particles N, whose sizes are distributed according to the 
desired N{a). It practice, however, it transpires that the CE represents a far from optimal 
framework for simulating polydisperse fiuids. The principal difficulty lies with the limited 
range of computationally accessible particle numbers which, in any simulation, is typically 
many orders of magnitude smaller than found in an experiment. The resulting finite-size 
effects are particularly pronounced in a CE simulation because the specific realization of the 
disorder N{a) is fixed. This suppresses large scale fiuctuations in p{a) and potentially leads 
to sampling deficiencies 0] . Additionally, the CE suffers other drawbacks familiar from sim- 
ulation studies of monodisperse fluids. For instance, relaxation times are extended because 
density fluctuations decay solely via diffusion; there is no direct access to information on 
chemical potentials; metastability and hysteresis hinder the study of phase transitions. 

Experience with the simulation of monodisperse fluids has shown that use of the grand 
canonical ensemble (GCE) is highly effective in circumventing many of the aforementioned 
problems associated with the CE |]5|-^]. As we shall show, its application in the context 
of polydisperse fluids retains many of the beneflts of the monodisperse case. Moreover, it 
provides the key to improved sampling of the density distribution p(o"). This is because 
within the GCE framework N{a) fluctuates as a whole, thereby capturing fluctuations in 
p(cr) on all simulation length scales. Notwithstanding these advantages, however, the GCE 
might appear, at flrst sight, unsuitable for the purpose of traversing a particular trajectory 
through the space of the density distribution p(o"). This is because p{a) ostensibly lies out- 
with the direct control of the simulator, its form being instead determined by the imposed 
chemical potential distribution fi{cr). Nevertheless, it turns out to be possible to tune fi{cr) 
within a histogram extrapolation scheme, in such a way as to realize a speciflc desired form 
of p(cr). We shall demonstrate that this dual use of the GCE and histogram extrapolation 
methods permits a chosen phase space path to be followed efficiently. 

The layout of our paper is as follows. In sec. |11 A| we formulate the statistical mechanics 
for a polydisperse fiuid within the grand canonical ensemble. We then describe (sec. [II B|) 
the combined GCE plus histogram extrapolation methodology for tracking a particular path 
through the space of p(o"). In section |ITI| we apply the method to the problem of obtaining 
the dilution line properties of three types of size-disperse hard sphere fiuids. The chemical 
potential distribution of these fiuids is determined as a function of volume fraction and the 
results compared with the predictions of two commonly used equations of state. Finally in 
section we discuss our findings and their implications. 



II. METHOD 
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A. Statistical Mechanics 



We consider a classical fluid of polydisperse particles confined to a volume V = L'^. The 
system is assumed to be thermodynamically open, so that the particle-number distribution 
N{a) is a statistical quantity. The associated grand canonical partition function takes the 
form: 

^v=T. JJ^Iil dr J da, exp {-^Hn ({r, a})) (2.1) 

with 

N 

Hn ({r, a}) = $ ({f, a}) - J] /i(a.) . (2.2) 

i=l 

Here N is the overall particle number, while (3 = {kBT)~^ and yu(o") are respectively the 
prescribed inverse temperature and chemical potential distribution, {r, a} denotes the con- 
figuration, i.e. the complete set (rl, ai), (r2, (72) • • • (rjv, cr^r) of particle position vectors and 
polydisperse attributes. The corresponding configurational energy $ ({r, a}) is assumed to 
reside in a sum of pairwise interactions 

N 

* {{r, cr}) = ^ 0(1 rl - f, |, a„ a,) , (2.3) 
i<j=i 

where is the pair potential. 

The instantaneous number distribution is defined by 

TV 

Ar(a) ^ ^ 5(a - a.) , (2.4) 

i=l 

with a the continuous polydispersity attribute and N = J N{a)da. We shall be concerned 
with the fluctuations in the associated density distribution 

p{a) ^ N{a)/V , (2.5) 

and the configurational energy density 

u = ^{{r,a})/V . 

The statistical behavior of these observables is completely described by their joint probability 
distribution 

VvW)M = {5(u- V-^^{{T,a}))\{5{p{a) - V~'N{a))^ , (2.6) 
or more explicitly 

PvW)M = ^T.wE\ daA exp {-PHn ({r, a})) (2.7) 



X 
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Integrating over the energy fluctuations yields the probabihty distribution function of the 
density distribution 



Pv[p{(y)] = J Pv[p{cr),u]du . (2.8) 
Our speciflc concern is with the average form of p{cr), given by 

p{a) = [ p{a)pv[p{cr)]dp{a) . (2.9) 



Given a prescribed chemical potential distribution /i(cr) and temperature j3, the form of 
p(cr) can be determined by simulation. Except in the ideal gas limit, however, no exact 
relationship between p{cr) and p{(t) will generally be available. Thus one cannot (from a 
'bare' GCE simulation), readily determine that /i(cr) corresponding to a particular desired 
target density distribution pt{cr). Subject to certain restrictions however, this can be achieved 
via use of histogram extrapolation. 

The key idea of histogram extrapolation is that a measured distribution [p(c"), m] ac- 
cumulated at one set of model parameters fi{a), P can be reweighted to yield estimates of the 
distribution appropriate to other parameters fi'{a),(5'. In its simplest form the reweighting 
is given by 

p'y[pia),u\p'ia),/3']=wpv[picr),u\pia),/3] (2.10) 
where the reweighting factor w takes the form 



N 



W 



exp Wp'{(^^) - /5/i(^.)] -V{P'-/3)u] (2.11) 



u=l 



By tuning the form of p'{(j) and the value of (3' within the reweighting scheme, it is possible 
to scan the space of p(cr), thereby "homing in" on the target density distribution. To this 
end it is expedient to deflne a cost function measuring the deviation of p(a \ p'{a),P') from 
the target form 

A{p\a),P') ^ I I p(a) - pt{a) I ^{a)da . (2.12) 

Here, for numerical convenience, we have incorporated a weight function 7(0") into our 
deflnition, the role of which (as described below) is to ensure that comparable contributions 
to the cost function arise from all sampled regions of the cr-domain. Within this framework, 
the task of determining those values of p\a),[3' that yield the target distribution pt(cr) 
reduces to that of functionally minimizing the cost function A with respect to p'{a),j3'. As 
we describe in the next section, this is achievable using standard algorithms. 



B. Implementation 

We have employed Monte Carlo (MC) simulation within the grand canonical ensemble 
to study the dilution line properties of systems of size-disperse hard spheres. This section 
details the principal aspects of our simulation and analysis procedure. 
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1. Program, data structure and acquisition 



Our simulated system comprises a variable number of hard spheres contained within 
a periodic box of volume V = L'^. The dimensionless polydisperse variable a was taken 
to be the particle dmmeter expressed in units of the mean diameter (see section pTI) . An 
upper bound was placed on the permitted range of diameters, and the simulation volume 
was partitioned into an array of cells each of linear dimension cTc, so that L = lac- 
This strategy aids efficient identification of particle interactions by ensuring that interacting 
particles occupy either the same cell or directly neighboring ones. 

The grand canonical ensemble Monte Carlo (GCMC) algorithm employed has a Metropo- 
lis form and invokes four types of operation: particle displacements, particle insertions, 
particle deletions and particle resizing; each is attempted with equal frequency. Specific 
to the polydisperse case is the resizing operation which entails attempting to change the 
diameter of a nominated particle by an amount drawn from a uniform random deviate con- 
strained to lie in some prescribed range. This range (maximum diameter step-size) is chosen 
to provide a suitable balance between efficient sampling and a satisfactory acceptance rate 
at the prevailing number density. As regards the remaining types of moves, these proceed 
in a manner similar to the monodisperse case except that for insertion attempts the 
new particle diameter is drawn with uniform probability from the range a e [0, Uc] . 

As primary input, the program takes the chemical potential distribution /i(o"), which is 
required for the accept/reject Monte Carlo lottery. This distribution is stored in the form 
of a histogram, constructed by dividing the truncated interval < a < ac into a prescribed 
number M of sub-intervals or "bins" [ill]. Thus all particles whose diameters fall within the 



scope of a given bin are associated with the same value of the chemical potential. 

The principal observable of interest is the joint probability distribution py[p{a),u] (c.f. 
Eq. p.8| ). Operationally it is impractical to construct the full form of this distribution in a 
simulation because to do so would entail forming a histogram over histograms-the memory 
storage costs of which would be prohibitive. The procedure adopted, therefore, was to 
sample the instantaneous density histogram p{a) (discretized in the same manner as /^(ct)), 
and append successive measurements of this quantity to a file |]T^ . The set of all such samples 
constitutes a list representing a sequential history of the individual data measurements [0 . 
This data list is post-processed by an analysis program which reads in each of the individual 
list entries and averages over all of them to construct a histogram approximation to the 
average density distribution p(o"). If desired, the analysis program additionally implements 
histogram reweighting of the data in order to enable extrapolation to neighboring values 
of the model parameters. This extrapolation is achieved by assigning a weight to each list 
entry of the form given by equation |2.11| . The complete set of weights permits construction 
of the reweighted histogram. 



2. Tracking a phase space path 

The strategy we have adopted is to traverse the phase space path of interest in a stepwise 
fashion, utilizing histogram extrapolation to proceed from one step to the next. In the 
general case of particles having a finite potential, the phase space path may involve changes 
in the temperature as well the form of p(o"). For simplicity of illustration, however, let us 
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presume that the path is isothermal (i.e (3 = constant) and that the form of the chemical 
potential distribution, /i'-^-' (a) say, is known at some arbitrary point along the path. 

The procedure is then as follows. From simulations at the known state point, data for 
Pv[p{cr),u] can be accumulated directly. Histogram reweighting is then applied to this data 
to extrapolate some distance along the path to a new point p^^\a), and to provide an 
estimate of the corresponding form of p^^\a). The latter quantity is then employed in a 
fresh simulation, the results of which are extrapolated to a point further along the path, and 
so on. By iterating this procedure p^*^(cr) — > p^^'^^\a), p^''\cr) — > p^'^'^^^a), one traces out 
the entire phase space path. 

The implementation of the extrapolation stage necessitates a prior choice for the step 
size, that is the difference between the measured p(cr) and the next target pt(cr). The 
magnitude of this difference should be chosen to be as large as possible, consistent with 
remaining within the range of reliable extrapolation. A good indicator that this is in fact 
the case is that the individual densities of the target distribution pt{cr) each overlap with 
the range of typical fluctuations appearing in the simulation data p{a) 

Once a suitable step size has been determined, the extrapolation procedure proceeds 
by minimizing (within the reweighting scheme) the cost function A introduced in eq. |2.12| . 
For all but the lowest densities, this task is complicated by the existence of strong coupling 
between the p variables, deriving from the fact that the number density for each a depends 
on the whole chemical potential distribution. Fortunately, efficient algorithms for performing 
multi-dimensional functional minimization are widely available and, at least for the cases 



we have considered, appear to operate effectively. The sole difficulty encountered was that, 
on occasion, the minimization failed to fully converge for a values in the wings of pticr)- 
To remedy this problem, a weight function 7(0") was incorporated in the cost function (c.f. 
eq. |2.12| ) , the purpose of which is to enhance the contribution to A from a values for which 
Pt{cr) is small. Good results were obtained by setting 7(0") oc [pt{cr)]^^ |jl5|. 

It should be stressed that the methodology set out above presumes the availability of 
a form of p(o") for some starting point on the phase space path of interest. This may be 
obtained straightforwardly if the path passes through a region of low density where reliable 
analytical estimates for p(o") can be employed. Otherwise the method must be bootstrapped 
by other means. One simple but effective approach for achieving this operates as follows. 
Starting from some initial guess for the desired p{a) (e.g. p°(cr) = ln[pt((T)]) a series of short 
simulations are carried out in which p{cr) is iterated according to 




P" 

where < 5 < 1 is a damping factor, the value of which may be tuned to optimize conver- 
gence. Although one could certainly envisage more sophisticated and efficient schemes, we 
have found this method to operate satisfactorily. 



III. DILUTION LINE STUDIES OF POLYDISPERSE HARD SPHERE FLUIDS 
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A. System and simulation details 



We have obtained the dilution hne properties (cf. sec. |) of size-disperse hard sphere |]TB 
fluids with diameter distribution N{a) assigned one of two forms: 

(i) Schulz distribution. 

(ii) Log-normal distribution. 

These distributions are conveniently expressed in terms of a normalized size function n{(T) = 
N{a)/N. For the Schulz, this takes the form 



I /z+1 



z+l 



z + l 



a 



a 



(3.1) 



where a is the average particle diameter and 2; is a parameter controlling the width of the 
distribution. For the log-normal distribution, the size function is given by 

aj2i,Hl + W^) \ 21ri(l + H'-^) 



with W the standard deviation in units of a. Note that both of these distribution are 
normalized, that is n{a)da = 1, and vanish as a — > 0, implying a natural lower limit to 
cr. By contrast, there is no finite upper limit and consequently, for simulation purposes, it 
was necessary to impose an upper bound (cut off) Uc (see also sec. |11 lj| ). 

We have studied three distinct size distributions — two of the Schulz form and one of the 
log-normal form. For the Schulz distribution, width-parameter values of 2; = 15 and z = 5 
were considered, with cutoff values cJc = 3 and cTc = 4 respectively. For the log-normal 
distribution, the single case W = 2.5 with Uc = 12 was studied. In each instance we set 
a = 1. Cell array sizes (cf. sec. P B 1|) of linear dimension / = 3, 4, 5 were used. In absolute 



dimensionless units (L = lac) these correspond to L = 9, 12, 15 respectively for the z = 15 
Schulz, to L = 12, 16, 20 for the z = 5 Schulz, and to L = 36, 48, 60 for the log-normal 
distribution. The histogram discretization parameter (cf. sec. [II B|) was set to M = 75 and 
M = 100 for the z = 15 and z = 5 Schulz distributions respectively, and to M = 120 for the 
log-normal distribution. 

The average density distribution can be expressed in terms of the normalized size function 
multiplied by the overall number density po = ^/V, i.e. 

p{a) = pon{cr) . (3.3) 

The procedure for tracking the dilution line utilized this overall density as a measure of 
the location along the line. The tracking procedure was initiated at a small value of po 
by approximating the chemical potential distribution according to the ideal gas relation 



/i(cr) = lnp(o") |T^. Histogram extrapolation was applied to the resulting simulation data 



in order to refine the initial estimate of p(cr). Thereafter the strategy described in sec. |11B 2 



was implemented to follow the dilution line to higher densities. At each step this entailed 
setting a target form for p(cr) corresponding to a value of po larger than that used at the 
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previous step. The cost function measuring the discrepancy between p{a) and the target was 
then minimized to yield an estimate for the appropriate On efficiency grounds, this 

minimization was performed in two stages; an initial approximation to was obtained 
from a one-dimensional minimization in which the activity distribution exp(/x((j)) was mul- 
tiplied by an overall factor. Thereafter the complete functional minimization was performed 
to yield a more accurate form for /u(cr). 

Although (within the specific context of the tracking procedure), po provides a convenient 
measure of the location on the dilution line, it conveys little information regarding the degree 
of packing within a polydisperse system. We shall therefore find it convenient to quote values 
for the overall volume fraction of the system given by 



It is this quantity (rather than po) which features in the presentation of our results, to which 
we now turn. 



Owing to the computational complexity of the simulations, we have obtained the com- 
plete dilution line for each fiuid only for the / = 3 cell array size. For the larger system 
sizes, only a few spot measurements were made along the dilution line, each of which was 
bootstrapped by appeal to the measured fi{a) obtained from the / = 3 systems. 

The dilution lines for the three fiuids were tracked to the highest computationally ac- 
cessible volume fraction, terminating only once the relaxation timescale became excessively 
large. The maximum volume fractions attained were t] = 0.445 and t] = 0.40 for the Schulz 
density distributions with z = 15 and z = 5 respectively, and r] = 0.33 for the log-normal 
case. Snapshot configurations for all three fiuids are shown in fig. |l| for rj values slightly 
below the maximum attained in each case. We mention in passing that, at least for the 
cases of the two Schulz distributions, the maximum volume fractions reached are somewhat 
larger than would be readily attainable for GCMC simulations of monodisperse hard spheres. 
The latter, of course, become highly inefficient at large rj because the acceptance rate for 
particle insertions falls rapidly as the free volume diminishes. Whilst the same is true in the 
polydisperse context for insertions of large particles, space can often be found for placing 
a small particle. This facilitates fiuctuations in the overall particle number A^, while resiz- 
ing operations (whose efficiency at large rj is typically greater than that of inserting large 
particles) ensure continued proper sampling of the density distribution. 

The measured forms for p(o") as a function of rj differ qualitatively between the Schulz and 
the log-normal density distributions, and accordingly we discuss them separately. Beginning 
with the Schulz case, figs. ^ and | show the measured p{a) and the corresponding /^(cr) for the 
two fiuids at a selection of volume fractions along their respective dilution lines. The range 
of a values shown is that for which the simulations delivered data of reasonable statistical 
quality. One notes that for both the z = 15 and z = 5 cases, the tails marking the large-cr 
vestiges of the distributions fall considerably short of the respective cutoffs a^- Indeed, we 
have verified that in the course of the simulations, no particles of diameters approaching the 
cutoff diameter occurred on the simulation timescale, implying that our data is unaffected by 




(3.4) 



B. Results 
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its imposition. As regards the forms of the chemical potential distributions, one sees that for 
small 7] they display a maximum near the peak in — behavior that is of course mandated 
at sufficiently low r] by the known properties of the ideal gas limit. For larger t], however, the 
peak is lost and /i(cr) increases monotonically. The increase of in the regime of large a 
indicates that the excess chemical potential (measuring the work associated with inserting a 
sphere) grows faster with a than the decrease in the ideal contribution associated with the 
decay of p(cr). 

Turning now to the log-normal case, fig. |^(a) demonstrates that p{cr) decays extremely 
slowly with increasing a. The peak in the distribution therefore occurs at much smaller a 
than for the two Schulz forms although the average diameter a is identical in all three cases. 
As a practical consequence of the slow decay (and notwithstanding the imposition of a very 
large value of dc), ranges of particle diameters extending up to the cutoff value were observed 
in the system. Indeed it was not feasible within the computational constraints to utilize a 
cutoff for which this did not occur, and hence truncation effects are always significant in 
this system — a point to which we shall return below. The measured forms of /^(cr) for the 
log- normal fiuid are shown in fig. ^(b). In contrast to the Schulz distributions, they display 
(for all accessible rj) a narrow maximum close to the peak in p{a) at small p. Thereafter with 
increasing a there is a slow fall to a broad minimum, whereafter p{(j) increases strongly. 

Having outlined the main qualitative features of the relationship between p((t) and p{a), 
it is instructive to perform a detailed comparison between our measurements and the pre- 
dictions of analytical equations of state (EOS) appearing in the literature. For hard spheres, 
two commonly used equations are that due to Boublik, Mansoori, Carnahan, Starling and 
Leland (BMCSL) [l^Jl^ based on the Carnahan-Starling equation for monodisperse hard 



spheres, and that due to Salacuse and Stell ||20|, based on the Percus-Yevick (PY) theory. 



Both are reproduced in appendix and express p{(j) in terms of an expansion to third 
order in a, with coefficients given in terms of the first three moments of p(cr). We have 
compared the predictions of the EOS for each of the three fiuids studied, with the finite-size 
simulation data at three values of r], namely a low, a moderate and a high value. We describe 
our findings for each fiuid in turn. 

The results for the Schulz distribution with z = 15 at the low volume fraction t] = 0.056 
are shown in fig. ^(a). At first sight there is good agreement between the L = 9 and L = 12 
simulation data and both the BMCSL and PY equations of state over the entire region of a. 
However, closer inspection reveals appreciable discrepancies between theory and simulation, 
not visible on the scale of p{a). These are apparent if one suppresses the dominant ideal 
gas contribution to expose the excess chemical potential, given by (see appendix ^ 

PcM) = /^(^) - ln[po'^(o-)] • (3.5) 

This quantity is plotted in the inset of fig. ^(a), from which one sees that compared to the 
simulation results, the EOS slightly underestimate p{a). 

Fig. ^(b) shows the results for the z = 15 Schulz distribution at the moderate volume 
fraction rj = 0.257. Again there is good agreement between the L = 9 and L = 12 simu- 
lation data suggesting that finite-size effects are insignificant. Here, however, discrepancies 
between the BMCSL and PY equations of state are larger than at the lower value of 1], being 
evident on the scale of the absolute chemical potential. One sees that both EOS significantly 
underestimate p{cr) for large a, although they agree quite well with one another. A simi- 
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larly high level of agreement between the data from the L = 9 and L = 12 system sizes is 
manifest at the higher volume fraction t] = 0.426, fig. ^(c). Here the PY equation of state is 
seen to fare somewhat better than the BMCSL equation although both underestimate fi{cr) 
substantially towards the upper end of the a range. 

A similar picture emerges for the Schulz distribution with z = 5 (fig. Again both 
equations of state underestimate fi{a), even at the lowest volume fraction, although the 
PY equation corresponds significantly more closely to the simulation results at high volume 
fractions than does the BMCSL equation. We could again discern no evidence for appreciable 
finite-size effects in the simulation results. 

In seeking to compare the results for the log-normal system with the EOS predictions it is 
essential to bear in mind the importance of truncation effects. The moments of a truncated 
log-normal distribution can differ dramatically from those of the full distribution even for 
large values of cXc. In order to facilitate a fair comparison with theory, the analytic form 
of /i((T) must therefore be calculated using the moments of the same truncated distribution 
as employed in the simulations. The results of performing this comparison are presented in 
fig. 1^. At low volume fraction, there is good agreement between the EOS predictions and 
the simulation results, while at high volume fraction the EOS underestimate the measured 
yu(cr) considerably, with the degree of discrepancy increasing towards the tail of the distri- 
bution. Once again we could discern no evidence of finite-size effects within the statistical 
uncertainties of our data. 

The results presented above indicate that the BMCSL and PY equations fail to provide a 
quantitatively accurate description of the chemical potential distribution, particularly when 
the volume fraction is large. It is instructive to examine the implications of this finding 
for calculations of densities, which depends very sensitively (indeed exponentially so) on 
the chemical potential. To this end we have studied the degree to which the form of /i((j) 
calculated via the BMCSL equation from some prescribed p(cr), actually yields this density 
distribution when input to a simulation. The results of performing this comparison are 
shown in figures §(a) - §(c) for each of the three fluids at a high volume fraction. In 
each instance, the solid line shows the input density distribution p{a) from which yu(cr) is 
calculated. The data points show the simulation results obtained using this form of /i(cr) 
for 3 different system sizes. As the figures clearly demonstrate, the measured form of p{a) 
deviate substantially from the prediction. 

Finally in this section, we examine the moment structure of the excess chemical poten- 
tial /iex(c") given by eq. P75|. Both the BMCSL and PY equations assume that He^i^a) is 



expressible in terms of a cubic polynomial in a. Bartlett has reached the same conclusion 
using geometrical arguments inspired by scaled particle theory We have investigated 
this proposal by fitting our data to the expression 

Aiex(o-) = - ln(l -rf) + aio + a2(y'^ + a^^a^ , (3.6) 

where the constant term is fixed by the requirement that in the limit a — > the probability 
of inserting a sphere is proportional to 1 — r/. We find that all our data are fitted very well 
by this expression; fig. ^(a) shows a typical fit for the case of the Schulz fluid with z = b 
a.t 1] = 0.377. Also shown, (fig. |^(b)) are the fit coefficients ai, 02,0:3 for various values of 
1], together with the predictions of the BMCSL equation of state. One sees that at high 
volume fraction, the BMCSL underestimates all coefficients, the relative discrepancy being 
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marginally larger for the ai coefficient than for the others. 



IV. DISCUSSION AND CONCLUSIONS 

In summary, we have presented a grand canonical simulation method for studying poly- 
disperse fluids. The method utilizes histogram extrapolation techniques to track efficiently 
an arbitrary path of interest through the space of the density distribution p(o"). We have 
applied it to the specific problem of obtaining the dilution line properties of size-disperse 
hard spheres. It should also prove useful in studying more general species of polydisperse 
fluids and their phase transitions, both in the bulk and confined geometries. We intend to 
report on such extensions in future communications. 

Previous simulation studies of polydisperse fluids have generally operated within a semi- 
grand canonical framework in which a fixed number of particles are studied either at 
constant pressure (see eg. ref. [^) or within a Gibbs Ensemble MC scheme |2^. In common 



with the present work, these studies utilized a fluctuating particle size distribution, realized 
by means of MC resizing moves controlled by a chemical potential distribution. In contrast 
to our approach, however, thermodynamic properties were studied as a function of the shape 
of the activity distribution exp(/i(cr)); no constraints were placed on the conjugate density 
distribution, which was consequently free to adopt whichever functional form minimized the 
free energy for the imposed /i(cr). In view of this, one might question the extent to which the 
simulation results reflect the actual situation in realistic systems. One situation in which 
the lack of a constrained density distribution might be relevant is the interesting issue of 
the influence of polydispersity on the freezing of hard spheres P3. This was investigated 



in ref. ||2^, using semi-grand canonical MC and Gibbs-Duhem integration. It would be of 
considerable interest to harness the approach described in the present work to investigate 
the effects on the freezing behaviour of varying the width of the density distribution whilst 
constraining its shape to some physically realistic form. 

In the present study, attention was focused on fluids having a high degree of polydisper- 
sity. The motivation for this choice was two-fold. First, models exhibiting a wide density 
distribution provide a suitably testing challenge against which to assess the effectiveness of 
our method. Second, there exist in the literature a number of interesting predictions con- 
cerning the role of depletion forces in highly polydisperse systems. For instance, it has been 
suggested by several authors that attractive depletion forces might engender novel phase 
transitions in polydisperse hard spheres P^^5|. Specifically, Sear has suggested that a fluid 
of hard spheres having a log-normal size distribution will be unstable with respect to crys- 
tallization of the large particles at a// finite volume fractions [^. By contrast, Cuesta 



has predicted that a sufficiently wide log-normal distribution will exhibit fluid-fluid phase 
separation at some finite density. It seems likely that for a truncated size distribution a 
phase transition will not occur for arbitrarily small volume fraction because of the absence 
of the largest particles which mediate the greatest depletion forces. Thus our simulation re- 
sults for the truncated distribution are able neither to conclusively conffim nor refute these 
predictions. It suffices to say that we observed no evidence of crystallization in the particular 
truncated log-normal fluid studied, up to a volume fraction of 77 = 0.33. Similarly, within 
the range of accessible volume fractions, no evidence for phase transitions was observed in 
either of the two Schulz fluids studied. 
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Turning finally to the comparison between our simulation results and the predictions 
of the polydisperse equations of state, we find that neither the BMCSL nor PY equations 
offer a quantitatively accurate description of the thermodynamics of hard spheres for large 
polydispersity and at high volume fraction [^. Both equations underestimate fi{a) at all 
fluid densities over the entire range of a, implying that they underestimate the Helmholtz 
free energy density / = Jq° fi{a\pQn{a))dadp'f) and hence overestimate the stability of the 
fluid. The magnitude of this overestimate becomes more pronounced the greater the volume 
fraction. Interestingly we find that in this regime the PY equation performs appreciably 
better than the BMCSL equation despite the fact (cf. appendix ^ that the latter derives 
from a monodisperse hard sphere equation of state which has been found to be superior to the 
PY approximation. It remains to be seen to what extent the overestimate of fluid stability 
impinges on the results of existing calculations of depletion-force induced phase separation 
based on the BMCSL and PY approximations l2^j2^ . In any case, our results should provide 
a useful testing ground for any future improvements to the existing polydisperse equations 
of state pO . 
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FIG. 2. Dilution line properties of hard spheres having the Schulz density distribution (z = 15, 
(Jc = 3.0) for system size L = 9. (a) Data points show the measured density distribution p{a) at 
a selection of values of volume fraction rj along the dilution line; dotted lines correspond to the 
target distribution pt{<T) (b) The corresponding chemical potential distribution //(cr). Statistical 
errors do not exceed the symbol sizes. 
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FIG. 3. Dilution line properties of hard spheres having the Schulz density distribution {z = 5, 
o"c = 4), for system size L = 12. (a) Data points show the measured density distribution p((t) at 
a selection of values of volume fraction rj along the dilution line; dotted lines correspond to the 
target distribution pt{(j). (b) The corresponding chemical potential distribution iJ,{a). 
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FIG. 4. Dilution line properties of hard spheres having the log-normal density distribution 

(W = 2.5, (Tc = 12) for system size L = 36. (a) Data points show the measured density distribution 
p((7) at a selection of values of volume fraction t] along the dilution line; dotted lines correspond 
to the target distribution pticr). The inset shows the region of small a. Lines are merely guides 
to the eye. (b) The corresponding chemical potential distribution //(cr). Statistical errors do not 
exceed the symbol sizes. 
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FIG. 5. Chemical potential distribution for the Schulz density distribution (z = 15, 

(Tc = 3), for system sizes L = 9 and L = 12. (a) r] = 0.056. Here the inset shows the excess 
chemical potential given by eq. |3.5| (b) rj = 0.257; (c) i] = 0.426. Also shown in each case are the 
predictions of the BMCSL and PY equations of state. Statistical errors do not exceed the symbol 
sizes. 
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FIG. 6. Chemical potential distribution /u((t) for the Schulz density distribution (z = 5, cJc = 4). 
(a) T] = 0.048, L = 12, 16, the inset shows the excess chemical potential given by eq. |3.5| ; (b) 
r] = 0.231, L = 12, 16; (c) rj = 0.36 L = 12, 16, 20. Also shown in each case are the predictions of 
the BMCSL and PY equations of state. Statistical errors do not exceed the symbol sizes. 
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FIG. 7. Chemical potential distribution fJ-{cr) for the log-normal density distribution (W = 2.5, 
(Tc = 12) for system sizes L = 36,48,60. (a) t] = 0.126, the inset shows the excess chemical 
potential given by eq. |3.5| ; (b) rj = 0.307. Also shown in both cases are the predictions of the 
BMCSL and PY equations of state. Statistical errors do not exceed the symbol sizes. 
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FIG. 8. The measured density distribution p{a) obtained using forms of /^(ct) predicted by 
the BMCSL equation of state at the stated volume fraction, (a) Schulz {z = 15), rj = 0.427, 
L = 9, 12, 15; (b) Schulz {z = 5), rj = 0.36, L = 12, 16, 20; (c) log-normal distribution, rj = 0.307, 
L = 36, 48, 60, the inset shows the same data on a log scale. In each case the density distribution 
from which /u((t) derives is shown as a solid line. Statistical errors do not exceed the symbol sizes. 
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FIG. 9. (a) Data points show the measured excess chemical potential distribution iJLe^{(y) for 
the Schulz density distribution (z = 5) at = 0.377, for system size L = 12. The solid lines shows 
a fit to the data of the form given by eq. |3.6| . (b) Values of the fit coefficients ai,a2,a3 for a 
selection of values of the volume fraction rj. The solid lines show the predictions of the BMCSL 
equation. Statistical errors do not exceed the symbol sizes. 



APPENDIX A: POLYDISPERSE HARD SPHERE EQUATIONS OF STATE 

The equations of state that we quote here express the chemical potential as a function 
the sphere diameter a. The expression due to Salacuse and Stell is a generalization of 
the Percus-Yevick result for monodisperse hard spheres: 
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/ipY(o-) = ln[po'^(cr)J - ln(l - C ) H 



l-Cs 

+ T-Ca)^"'" + (Al) 



The BMCSL equation of state [|T8| , p!9| , on the other hand, generahzes the Carnahan-Starhng 
expression, which for monodisperse hard spheres is more accurate than the Percus-Yevick 
result. It is given by 

/iBMCSL(o^) = HPon{a)] + (Sa^^VCl - 2a\|/C|) ln(l - C') 
, cr3(Co- CI/CI) + 3cr2Ci + 3aC2 



1-C; 



3 



, ^^(3CiC2- CI/CI) + 3a^Cl/C3 , 2aXI .,,x 

where (n = T^Po^n/Q and m„ is the n-th moment of p(cr); note that = V: the volume 
fraction of hard spheres. 

Both equations of state include the ideal contribution to the chemical potential, 
/^idoai(o") = ln[po'^(<7)]- For general temperature — remember that we set /3 = 1 for our 
hard sphere system — this contribution would read /?/iideai(c) = Inp(cr). Because there has 
been some discussion in the literature regarding how the ideal chemical potential for poly- 
disperse systems should be assigned (see e.g. |2^), it may be helpful to note that, within 
the grand canonical framework, /iideai(<7) is unambiguously defined via the grand canonical 



partition function. Indeed, for an ideal system (Ti^r = 0), one readily finds from eq. |2.1| that 
p(cr) = exp(/3p(cr)), giving /3pideai(c") = lnp(o") as before. One may worry about dimensions 
in the argument of the log here; but the definition in eq. already implies that dimension- 
less units are used for lengths and for cr. Were this not the case, the integrals over fi and cTj 
would need to be normahzed by a unit volume vq and a unit value ctq for the polydisperse 
attribute, and the ideal chemical potential would read /5pideai(c") = ln[fo0'op(o")]) with the 
argument of the log now manifestly dimensionless. The normalization constant focxo could 
in fact be made dependent on cr; this would just give a a-dependent shift in the zero of the 
chemical potential scale. 
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